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[ This paper describes a new computer code which simulates a self-gravitating, 

relativistic perfect fluid in spherical symmetry. The fluid model uses an equation 
of state P = (r — l)p, where P and p are the fluid pressure and total energy density, 
respectively, and T is a constant satisfying 1 < T < 2. The code has been optimized 
for ultrarelativistic fluid flows, that is, for flows with Lorentz factors much larger than 
unity. This optimization involves a novel definition of the fluid variables, the use of 
a modern high-resolution shock-capturing scheme, and care in reconstruction of the 
primitive fluid variables — such as pressure and velocity — from the conserved quantities 
which are actually evolved by the code. 

Our new code was specifically developed to study the critical gravitational collapse 
of perfect fluids. Critical collapse has become an interesting subfield in general 
relativity since its initial discovery in the massless Klein-Gordon system jp, and the 
perfect fluid model has played an important role in advancing our understanding of 
the critical phenomena which arise at the threshold of black hole formation. (For an 
excellent introduction to critical phenomena, see the review by Gundlach Q.) While 
the critical solutions for perfect fluids in spherical symmetry have been the subject of 
recent study J3|, ^, |^, |(| [?J pi [| [n], , the precise nature of the critical solutions for 
r > 1.89 was not previously known, and thus one of the chief goals of our investigation 
was a thorough analysis of this regime. In the remainder of this paper we describe the 
equations of motion which are solved, and the numerical techniques which we use to 
solve them. A companion paper describes in detail the results we have generated 
with the code. 
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2. Geometry and fluid model 

The Einstein equations couple the spacetime geometry, encoded in the Einstein tensor, 
Gab, to the stress-energy tensor, T a b, associated with the matter content of the 
spacetime: 

Gab = 87rT afc , (1) 

(here and throughout, we use units in which the speed of light, and Newton's 
gravitation constant are unity: c = 1 and G = 1, and Latin indices a, b, c, • •• take 
on the spacetime values 0, 1, 2, 3.) A fluid is a continuum model for a large number 
of particles that uses macroscopic properties of a thermodynamic system, such as 
internal energy and pressure, as fundamental dynamical variables. A perfect fluid has 
no shear stresses or dissipative forces, and has a stress-energy tensor 

Tab = (p+ P)u a u b + Pg a b, (2) 

where p is the energy density, P is the pressure, u a is the fluid's four-velocity, and 
g a b is the spacetime metric. The energy density p contains all contributions to the 
total energy, which for a perfect fluid include the rest mass energy density, p a , and 
the internal energy density 

P = Po + Pot, (3) 

where e is the specific internal energy. The fluid number density, n, is related to p a 
via 

p = ran, (4) 

where m is the rest mass of a single fluid particle. The basic equations of motion for 
the fluid can be derived from local conservation of (a) the energy-momentum 

V a T ab = 0, (5) 

and (b) the particle number 

V a (nu a ) = 0, (6) 

where V a is the (covariant) derivative operator compatible with g a b- To these 
conservation laws one must adjoin an equation of state, P = P{p , e), which, further, 
must be consistent with the first law of thermodynamics. 

2.1. Equation of state 

The equation of state (EOS) closes the fluid equations by providing a relationship 
between the pressure and (in our case) the rest energy density and internal energy. 
The nature of this relationship provides much of the physics for a given system. As 
mentioned in the introduction, our primary motivation for exploring ultrarelativistic 
fluid dynamics is to study perfect fluid critical solutions. We expect these solutions 
to be scale invariant (self-similar), and we therefore choose an EOS compatible with 
this symmetry. The EOS 

p = (r - i) P , (7) 

where T is a constant, is the only EOS of the form P = P(p) which is compatible with 
self-similarity 13, O, 15 , and is notable for the fact that it results in a sound speed, 



c s , which is independent of density: 



Ultrarelativistic fluid dynamics 



3 



One can argue that this EOS is particularly appropriate for ultrarelativistic fluids, 
and hence we will refer to ([?]) as the ultrarelativistic equation of state. We note that 
the EOS for a "radiation fluid" corresponds to T = |, while T = 1 gives a pressureless 
fluid (dust). We do not consider the case of dust collapse here; hence, in what follows, 

i < r < 2. 

Another important fluid model is the ideal gas with the equation of state 

P = (T-l) Po e. (9) 

In the ultrarelativistic limit, the kinetic energy of the constituent particles of the 
fluid (or internal energy of the fluid in a thermodynamic context) is much larger than 
the mass energy, p e 3> p , giving p w p Q e. Thus, one can interpret the EOS (g) 



as the ultrarelativistic limit of the ideal-gas EOS. As discussed in 12 , the ideal- 
gas EOS, in the ultrarelativistic limit, becomes, in a limiting sense, scale invariant. 
As the critical solutions reside in this ultrarelativistic limit, the critical solutions for 
fluids with the ideal-gas EOS are reasonably expected to be identical to the critical 
solutions computed using (0). For this reason we hereafter limit our attention to the 
ultrarelativistic equation of state. 

2.2. Geometric equations of motion 

We use the ADM 3+1 formalism (specialized to spherical symmetry) to integrate the 
Einstein equations, and choose polar-areal coordinates for simplicity of the equations 
of motion and for singularity avoidance. Specifically, adopting a polar-spherical 
coordinate system (t, r, 9, <j)) 1 we write the spacetime metric as 

ds 2 = -a(r, t) 2 dt 2 + a(r, t) 2 dr 2 + r 2 {dO 2 + sin 2 d<t> 2 ) , (10) 

wherein the radial coordinate, r, directly measures proper surface area. In analogy 
with the usual Schwarzschild form of the static spherically symmetric metric, it is also 
useful to define the mass aspect function 

m(r,t) = lfl-iV (11) 



2 V a 2 . 

The fluid's coordinate velocity, v, and the associated Lorentz gamma function, W, are 
defined by 



au 



v(r,t) = — r, W(r,t) = au\ (12) 
cm 

Since the fluid four- velocity is a unit-length, time- like vector {u a u a = —1), we then 
have the usual relation between W and v: 

W 2 = - J-j. (13) 
1 — ir 

We now introduce two conservation variables 
r(r,t) = (p + P)W 2 - P 

S{r,t) = {p + P)W 2 v, (14) 
so named because they allow the fluid equations of motion to be written in conservation 



form (albeit with the addition of a source term), as discussed in detail in section 3.1, In 



contrast to the conservation variables, we refer to the quantities P and v as primitive 



Ultrarelativistic fluid dynamics 



4 



variables. With the above definitions, the non-zero components of the stress-energy 
tensor are given by 

T\ = -t T\ = Sv + P 

m _«c T\=T\ = P. (15) 



a 



<t> 



A sufficient set of Einstein equations for the geometric variables a and a are 
given by (a) the non-trivial component of the momentum constraint (the notation 
d x f denotes partial differentiation, i.e. d x f = df/dx) 

d t a = ~ATrraa 2 S, (16) 

and by (b) the polar slicing condition, which follows from the demand that metric 
have the form (10) for all t: 



d r Qna) = a 2 4irr (Sv + P) 
An additional equation for a(r, i), 
d r a = a 3 ( 47rrr 



rn 



rn 



(17) 
(18) 



follows from the Hamiltonian constraint. 



2. 3. Fluid equations of motion 

Given the ultrarelativistic EOS (]?]), the time evolution of our perfect fluid is completely 
determined by V a T ab = 0. The derivation of the equations of motion — which can can 
naturally be written in conservation form — is a straightforward piece of analysis, and 
will not be given in detail here. Instead, we will simply quote the results, and for 
convenience in discussing the numerical method of solution, we adopt a "state vector" 
notation. We thus define two-component vectors q and w, which are the conservation 
and primitive variables, respectively 

(19) 



(20) 
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' p ' 


q 


s 


w = 


V 



We then define a "flux vector," f , and a "source vector" ifi 



f = 



S 

Sv + P 



</> = 





s 



These variables have been introduced with a hat (*) to distinguish them from the new 
variables defined in section 2.4, which are subsequently used in the actual numerical 
solution algorithm. Further, to expedite the discretization of the equations of motion, 
we decompose the source term, S, into two pieces, as follows: 

(21) 
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ar 



where 



r ) [8iraarP 



aaP- 



(22) 



e = (Sv- 

\ r* J r*- 

We note that in spherically symmetric Minkowski spacetime we have = and 
S = 2P/r. With the above definitions, we can now write the fluid equations of motion 
in the conservation form 
1 



;0r 



(23) 
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where 

X = - (24) 
a 

is a purely geometric quantity. 

Written in the above form, the fluid equations of motion ( p^ ) contain a mixture 
of conservation and primitive variables, and thus it is necessary to transform between 
both sets of variables at each step in the integration procedure. The primitive 
variables w can be expressed in terms of the conservation variables q by inverting 
the definitions ( |l4| ) of the conservation variables: 

P = -20t + [4(3 2 t 2 + (T - 1)(t 2 - S 2 )] * (25) 
« = (26) 



where the non-negative constant [3 is defined by 

P = ~ (2 - T) . (27) 



The pressure equation ( |2q ) comes from the solution of a quadratic with a specific 
root chosen to yield a physical (non-negative) pressure. This demand (P > 0) further 
requires that r > \S\. A second physical requirement is that v be bounded by the 
speed of light, \v\ < 1, and from ( |26| ) this will clearly be automatically satisfied when 
t > \S\. These physical restrictions on the primitive variables can sometimes be 
violated in numerical solutions of the fluid equations, and we discuss some numerical 



techniques aimed at ameliorating such difficulties in sections 4.3 and 4.4. Finally, we 
note that the above transformation from q to w is particularly simple in that it can 
be expressed algebraically. The corresponding transformations for the gamma-law gas 
EOS (Q) involves a transcendental equation which, in a numerical implementation, 
must be solved iteratively at each grid point. 



2-4- New conservative fluid variables 

Using the conservation variables q defined above, and the numerical method described 
in sections || and ^, we developed a preliminary code to solve the relativistic fluid 
equations. We then tested this code by considering evolutions in Minkowski spacctime 
using slab and spherical symmetry. The tests in slab symmetry were completely 
satisfactory, modulo the convergence limitations of the numerical scheme. However, 
in spherical symmetry, we found that our method frequently failed for "stiffer" fluids 
(r > 1.9), most notably in "evacuation regions" where p — ► 0. Additionally, the fluid 
in such regions often became extremely relativistic, and the combination of p — > and 
1 17 j — > 1 proved particularly difficult to simulate. These problems that we encountered 
in spherical symmetry led us to seek a new set of conservation variables, and to 
motivate this change of variables, first consider the evolutions shown in figure [l]. Here 
we begin with a time-symmetric, spherical shell of fluid, which has a Gaussian energy 
density profile. Due to the time-symmetry, as the evolution unfolds, the shell naturally 
splits into two sub-shells — one in-going and one out-going — and as the sub-shells 
separate, a new evacuation region forms in the region where the fluid was originally 
concentrated. Examination of the conservation variable profiles reveals that l^l ~ r, 
and this observation suggests that we adopt new variables 

<I> = t-S, n = r + S, (28) 
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Figure 1. These plots show various fluid quantities at four different instances 
(equally spaced in time) in a flat spacetime, slab-symmetric evolution with 
r = f.9. The initial configuration is a time-symmetric Gaussian pulse. The 
top frames show the evolution of the original conservation variables, r and S. 
As the evolution proceeds, the pulse separates into left and right-moving halves, 
and a vacuum region (r — » 0) develops between the two sub-pulses. The bottom 
frames show the evolution of the new conservation variables, fl and <!>, which are 
specifically defined so as to avoid the formation of such vacuum regions. The 
correspondence of the new variables to left and right moving "waves" is also 
evident. Note that the plots of t, fl and $ have the same vertical scale, while the 
vertical scale for S is shown separately. The horizontal (radial) scale is the same 
for all of the plots. 



which loosely represent the in-going ($) and out-going (II) parts of the solution. Thus 
our new state vector of conservation variables is 

n 



(29) 



Not surprisingly, the numerical difficulties in evacuation regions are not completely 
cured with this change of variables; however, the new variables q provide a significant 
improvement over q in evolutions of spherically symmetric fluids with T > 1.9. 

The equations of motion for the new variables q can be readily found by adding 
and subtracting the two components of (p3|), giving 



1 



d t q+-d r [r 2 Xf] = t/>, 



(30) 
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where the flux and source terms are now given by 
i(II-$)(l + u)+P 



2 1 

!(n-$)(i-u)-p 











- £ 



(31) 



The transformation from conservative to primitive variables can be found by simply 
changing variables in and ( p6|) 



p = -/3(n + $) + 1 (n + $) 2 + (r - i)n$ 
n- $ 



(32) 
(33) 



n + $ + 2P 

We note that, given r > the new variables q are strictly positive: II > 0, $ > 



2.5. The perfect fluid as a scalar field 

There is a well-known relation between an irrotational, stiff (r = 2) perfect fluid and a 
massless Klein-Gordon scalar field. In this section we discuss the relationship between 
scalar fields and perfect fluids for < T < 2. The perfect fluid equations of motion 



VaT ab = 0, 

can be written in terms of p, P, and u a as 
u a V a + (p + P)V a u a = 0, 



(p + P)u a V a u b 



■9 



ab 



u a u b 



) V a P - 0. 



(34) 

(35) 
(36) 



If we assume the ultrarelativistic equation of state, P = (T—l)p, then these equations 
become 



v ( P 1/r V) 

u a V a u b 



0. 

(r-i) 
r 



(g ab + u a u b ) V a lnp = 0. 



(37) 
(38) 



We seek a specific combination of p and u a that allows the fluid equations to be written 
in terms of a single variable, and therefore introduce the ansatz 

(39) 



w 



where p is a constant that will be determined below. From elementary contractions 
we can express both p and u a in terms of w a 



p = {-w a w a )^ 



V 



-WbW 



w 



(40) 
(41) 



However, it remains to see if p can be chosen such that w a will satisfy the fluid 
equations of motion. We substitute expressions (|i(]) and (|d|) into the momentum 
equation pa), and find that this equation is satisfied provided that 



r- 1 



v = 



and 



[aWb] 



= 0. 



(42) 
(43) 
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This latter condition allows one to write w a as the gradient of a scalar field 

w a = V Q ^. (44) 

The equation of motion for <p is obtained from ( |35"| ) 

VJ(-V^V»' / VV]=0, (45) 

where 

v = W^ry (46) 

The condition (^), V[ a w b ] = 0, reduces to the requirement that the fluid be 
irrotational 

V [a u b] = 0. (47) 

Thus, the fluid equations for an ultrarelativistic, irrotational fluid can be written in 
terms of a nonlinear equation for a scalar field, if. For the stiff fluid (r = 2), we find 
that the equation of motion for (p becomes the massless Klein-Gordon equation 

V Q V<V = 0. (r = 2) (48) 

One typically places physically motivated conditions on the fluid variables, such 
as p > and u a u a = — 1. Solutions of the Klein-Gordon equation, however, have 
time-like, null, and space-like gradients (V a ip). With the usual physical constraints 
on the fluid, then only a subset of possible Klein-Gordon solutions can be interpreted 
as T = 2 perfect fluids, namely those with Vq^V'V < 0. 



3. Numerical methods for fluid equations 

An important consideration for numerical solutions of compressible fluid flow is 
how the numerical method will respond to the presence or formation of shocks, i.e. 
discontinuities in the fluid variables. These discontinuities often cause the dramatic 
failure of naive finite difference schemes, and as shocks form generically from smooth 
initial data, many special techniques have been developed for the numerical solution 
of fluid equations. One approach is to introduce an artificial viscosity that adds extra 
dissipation in the vicinity of a shock, spreading the would-be-discontinuity over a 
few grid points. This technique has been widely used, and has the advantages of 
simplicity of implementation and computational efficiency. However, Norman and 
Winkler jl6| investigated the use of artificial viscosity in relativistic flows, and showed 
that an explicit numerical scheme treats the artificial viscosity term inconsistently in 
relativistic fluid dynamics, leading to large numerical errors in the ultrarelativistic 
limit, W ^> 1. A second approach to solving the fluid equations with shocks comes 
from methods developed specifically for conservation laws. These methods, usually 
variations or extensions of Godunov's original idea [jlTj to use piece- wise solution of the 
Riemann problem, have proven to be very reliable and robust. LeVeque [l9| has 
written excellent introductions to conservative methods, and our presentation here is 
in the spirit of his work. Furthermore, the application of these methods to problems in 
relativistic astrophysics has been recently reviewed by Ibahez and Marti |2C| . However, 
for the sake of completeness, we first briefly define and discuss conservation laws, and 
outline a general approach for their solution. We then discuss a linear Riemann solver 
and a cell reconstruction method that results in a scheme which, for smooth flows, is 
second order accurate in the mesh spacing. 
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3.1. Conservation methods 

Conservation laws greatly simplify the mathematical description of physical systems by 
focusing on quantities Q — where Q may be a state vector with multiple components — 
that do not change with time 



dt dQ = 0. (49) 



v 



In this section we discuss the derivation of numerical schemes for this specific and 
important case where J dQ is conserved on the computational domain. Our discussion 



will be general, and not specifically tailored for the fluid PDEs derived in section 2.4, 
but for simplicity we restrict the discussion to one dimensional (in space) systems. 

While conservation laws are often written in differential form (e.g. V a T ab = 0) it 
is useful to first consider an integral formulation, which is often the more fundamental 
expression. Consider an arbitrary volume or cell, Ci, with a domain [xxj^a]- The 
quantity of Q within d is denoted Q;, and we define a density function q such that 

px 2 

Qi= dxq. (50) 

J X\ 

The change of Qi with time can be calculated from the flux, f(q), of q through the 
cell boundaries. This consideration thus yields our conservation law: 

d 



dt 



dxq(x,t) = f(q(x u t))-f(q(x 2 ,t)). (51) 

IXl 

The conservation law can be written in integral form by integrating ( pl|) from an initial 
time, t\, to a final time, t 2 , 

dxq(x,t 2 ) = 

Xl 

f 2 dxq(x,t 1 ) + f ' dtf{q(x u t))- [ * dtf{q(x 2 ,t)) (52) 

J XI Jti Jti 

and the differential form follows from further manipulation if we assume that q is 
differcntiable: 

d t q + d x t(q) = 0. (53) 

It should be emphasized that the integral formulation should be viewed as the primary 
mathematical form for a conservation principle, because it is not dependent on an 
assumption of differentiability. For example, at a shock front in a fluid system, q is 
not differentiable, and the differential form of the conservation law fails, while the 
integral formulation is still satisfied. Discretizations of conservation equations via 
finite differences rely on the differential form, and artificial viscosity must be added 
near shock fronts, forcing q to be differentiable. An alternate strategy is to develop 
numerical algorithms based directly on the integral formulation of the conservation 
laws. The Godunov method and its extensions are examples of this latter approach, 
and are the topic of the next section. 
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3.2. Godunov's Method 

Numerical algorithms for conservation laws are developed by discretizing the equations 
in their fundamental integral form. These methods derive from a control volume 
discretization, whereby the domain is divided into computational cells, Ci, now defined 
to span the interval [as — Ax/2, x + Ax/2] = [a;j_i/2, 2:1+1/2]; where Ax is the (local) 
spatial discretization scale. Following the derivation of the integral conservation 
law ( [52|) for the computation cell Ci, we introduce the averaged quantities, q": 
1 rXi+1/2 

q? = — / dxq(x,t n ), (54) 

with t n = nAt, where At is the temporal discretization scale. We then obtain the 
discrete form of the conservation law ( |52| ) 

qr +1 = qr-^(F l+ i/ 2 -F. t _ 1/2 ), (55) 

where the "numerical flux" is defined by 

Fi+i/2 = -^J dtf{q(x i+1/2 ,t). (56) 

At first blush, a numerical method based on a discretization of the integral 
conservation law does not appear promising: the flux integral ( |56| ) does not appear 
readily solvable, and it generally is not. However, in his seminal work, Godunov []l7f 
devised a technique to approximately evaluate the flux integral by replacing the 
function q(x,i„) with q[x,t n ), where q{x,t n ) is a piece-wise constant function. In 
this approach, the individual cells ("control volumes") arc treated as a sequence of 
"shock tubes", and a separate Riemann initial value problem is solved at each cell 
interface. Provided that the waves from neighboring cells do not interact — a proviso 
which gives a Courant-type condition on the time-step — each Riemann problem can 
be solved exactly to yield the local solution q(cc, t) (for t > t n ) for each "shock tube." 
Furthermore, since the solution of each of the local Riemann problems is self-similar, 
q(:Ei+i/2, £) is a constant in time, and the evaluation of the integral ( |56| ) becomes 
trivial. This then allows one to find explicit expressions for the cell averages at the 
advanced time, q n+1 , via (p35|). In summary, the Godunov method proceeds as follows: 
(a) From the average q™, one "reconstructs" a piece-wise constant function q(x, i„) to 
approximate the solution in Cf, (b) the Riemann problem is solved at the interfaces 
between cells, giving the solution q(x,t) for t n < t < t n+ i, (c) the solution q(x,t n +i) 
is averaged over the cell Ci to obtain the average at the advanced time, q™ +1 . We 
note that methods for solving the Riemann problem exactly for relativistic fluids have 
been given by Smoller and Temple [Q for the ultrarelativistic EOS, and by Marti and 
Miiller [|| for the ideal-gas EOS. 

Godunov's method has many nice properties: in particular, it is conservative and 
allows for the stable evolution of strong shocks. However, the original scheme does 
have some shortcomings: convergence is only first order, and the exact solution of the 
Riemann problem may be computationally expensive, especially for relativistic fluids. 
The convergence of the scheme can be improved by providing a more sophisticated 
reconstruction q(x,t n ), giving what are known as high-resolution shock- capturing 



methods. One such procedure is described in section 3.3, with details concerning 



the scheme's convergence given in section 4.7. In order to address the issue of 



computational efficiency, approximate Riemann solvers have been developed that 
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relate the problem-at-hand to a simpler system, for which the Riemann problem 
is easier to solve. Several approximate Riemann solvers have been developed for 
classical fluid dynamics, and many of these approximate methods have been extended 
to relativistic fluid systems. These include relativistic two-shock solvers [E3l E4J, a 



relativistic HLLE solver 1251] , and, as discussed in section 3.4, various linearized solvers 



3.3. Cell reconstruction 

Godunov-type numerical methods are based on solutions of the Riemann initial value 
problem at the interfaces between cells. As discussed above, during an update step one 
introduces functions q(x,i) — defined piece-wise on the intervals [£j_i/2j ^i+1/2] — to 
approximate the solution in the control volumes Ci. These functions are created from 
the cell averages q", and hence are called reconstructions. Consider the cell interface 
at x i+ i/ 2 : the state of the fluid immediately to the right (left) is q[ +1 / 2 (^+1/2)- The 
simplest reconstruction is to assume that q is piece-wise constant 

qf+1/2 = qi> qi+1/2 = qi+i, ( 57 ) 

as used in the original Godunov method and, as already discussed, this reconstruction 
results in a numerical scheme in which the spatial derivatives (and hence the overall 
scheme) have first order accuracy. The convergence can be improved by using a higher- 
order reconstruction for q, but care must be exercised so that the reconstruction does 
not induce spurious oscillations near discontinuities (see figure ||) . 

We have chosen to use a piece-wise linear reconstruction for q, which formally 
results in a scheme with second order con vergence. (The convergence properties are 



discussed in greater detail in section 4.7.) The q are reconstructed using the total 



variation diminishing (TVD) minmod limiter introduced by van Leer |26|. The van 
Leer limiter forces q to be monotonic near discontinuities, and this reduces the (local) 
accuracy of the scheme to first order. The first step of the reconstruction algorithm 
involves the computation of the slope (derivative of the dynamical variable) centered 
at the cell boundaries 

Si+1/2 = ■ (58) 

r i+ i - r, 

A "limited slope" , cr, , is then calculated via 

o-i =mhmod(s j _ 1/ 2,s j+ i/ 2 ), (59) 
where the minmod limiter is defined by 

!a if |a| < |6| and ab > 
b if |6| < |a| and ab > (60) 
if ab < 0. 

Using the limited slopes, we evaluate q at the cell interfaces as follows: 

qf+i/2 = q* + o-j(n+i/2 - n) (6i) 

and 

qI+i/2 = Qi+i + fi+ifa+i/a - r i+i)- (62) 

Finally, if we are unable to calculate physical values for w and w r (a situation which 
can and does occur owing to the finite-precision nature of our computations) we revert 
to a piece- wise constant reconstruction for q e and q r . 
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Figure 2. The three frames of this plot show different ways a discretized function 
can be reconstructed in a control- volume numerical method. The solid line 
represents a continuum (or "analytic" ) function and the solid hexagons represent 
discrete, approximate values of the function defined at grid points. Frame (a) 
represents the piece-wise constant reconstruction. Frame (b) shows a naive 
piece-wise linear reconstruction of each cell using s i+1 /2- This reconstruction 
oscillates near discontinuities in the function — such oscillations can easily lead 
to instabilities. Frame (c) shows a piece-wise linear reconstruction performed 
with the minmod limiter as described in the text. This reconstruction produces 
a discrete representation of the dynamical variable which remains well-behaved 
near discontinuities. 



3.4- The Roe linearized solver 

Perhaps the most popular approximate Riemann solver is the linearized solver 
introduced by Roe p7fl . This solver (and subsequent variants) has been used in a 
variety of applications involving general relativistic fluids |2^, ||^, |3(], 31, [32| [33|, and 
has proven to be robust and efficient. (The efficiency comparison is relative to solving 
either the exact Riemann problem for relativistic fluids, or a nonlinear approximation, 
such as the two-shock solver.) As the name suggests, the linearized solver approximates 
the full nonlinear problem by replacing the nonlinear equations by linear systems 
defined at each cell interface. The associated linear Riemann problems can then 
be solved exactly and cheaply, and the resulting solutions can be pieced together 
to produce an approximation to the solution of the original, nonlinear equations. 
Thus, in order to understand the Roc scheme, it is instructive to first consider linear 
conservation laws. 

The linear, scalar advection equation 

d t q + X8 x q = 0, (63) 

has the well-known solution q(x,t) = q(x — At, 0), where A is a constant and q(x,0) 
specifies the initial state. This scalar solution can be extended to linear systems of 
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conservation equations 

d t q + Ad x q = 0, (64) 

where A, an M x M constant matrix, is, by assumption, diagonalizable, with real 
eigenvalues, A M . (Greek indices take the values 1, . . . , M.) Let R be the matrix of 
right eigenvectors, r M , of A: 

R= [n|...|r M ], (65) 

and let A be the diagonal matrix: 

A = diag[Ai, . . . , Am]- (66) 

We then have 

A = RAR-\ (67) 

and the solution of the system may be obtained by introducing "characteristic 
variables", v: 

v = iT x q. (68) 

Using characteristic variables, the equations (|6~i| ) decouple into a set of scalar advection 
equations 

d t v + Ad x v = 0, (69) 

which can be immediately solved via: 

v M (x,i) =v M (x- A M t,0). (70) 

Given v(x, t), the transformation q = i?v then produces the solution of d64| ) in terms 
of the original variables, q. 

Turning now to the nonlinear case, the key idea is to first write the nonlinear 
system in quasilinear form 

d t q + A(q)d x q = 0. (71) 

Here, A is an M x M matrix which is now a function of q. Roe gives three specific 
criteria for the construction of A: 

(i) ^(q^q r )(q r -q*)=f(q r )-f(q*); 

(ii) «4(q , q r ) is diagonalizable with real eigenvalues; 

(iii) -4(q f , q r ) — > f (q) smoothly as q^, q r — > q. 

The latter two criteria can generally be satisfied by letting „4 be the Jacobian matrix 
evaluated using the arithmetic average of the conservation variables at the interface: 

0f(q) 



A = 



dq 

where 



(72) 

q=q i+ l/2 



q«+i/2 = 2 (qf+i/2 + qI+1/2) • (73) 



While this construction does not generally satisfy the first criterion, ( |72[ ) is often used 
in relativistic fluid dynamics (see for example J28|, |3(], on the basis of its relative 
simplicity, and we also adopt this approach. On the other hand, other authors |^9| have 
constructed a linearized Riemann solver for relativistic fluids with true Roe averaging, 
and we therefore refer to our scheme as a "quasi-Roe" method. 
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Having denned a specific linearization, the scheme proceeds by evaluation of 
_4(q i+1 / 2 ) — which is now viewed as a matrix with (piecewise) constant coefficients — 
followed by the solution of the Riemann problem for the resulting linear system. 
Carrying through an analysis not given here (see e.g. |18[), the Roe flux can be defined 

as 



_ 1 

-1/2 - g 



f(qm /2 ) + f(qW)-El-VI A ^ 



(74) 



-1/2) T HHi+1/2/ - 

M 

where, again, A M and r M are the eigenvalues and (right) eigenvectors, respectively, of 
■4(q i+ i/ 2 ). The quantities are defined in terms of the the jumps in the fluid 

variables across the interface 

q[+i/ 2 - qf+1/2 = Auvrv (75) 

For completeness, we give explicit expressions for the eigenvectors and eigenvalues of 



the ultrarelativistic fluid system (pOD in Appendix A 



Finally, it is important to remember that approximate Riemann solvers produce 
approximate solutions, which, under certain conditions, may diverge from the physical 
solutions. For example, concentrating on the Roe solver, Quirk p4| has recently 
reviewed several "subtle flaws" in approximate solvers. Fortunately, the approximate 
solvers often fail in different ways, and where one solver produces an unphysical 
solution, another solver may give the physical solution. Thus, it may be necessary 
to investigate a particular problem with multiple approximate Riemann solvers. 
Therefore, we have also implemented Marquina's solver |35]], an alternative linear 
solver that has also found application in relativistic fluid studies |3(| Q , as an option 
in our code. In addition to using the quasi-Roe and Marquina solvers to investigate 
the critical collapse of perfect fluids, we also implemented the HLLE solver in an 
independent code. We found that the quasi-Roe solver gave accurate solutions, and 
provided the best combination of resolution and efficiency for the critical collapse 
problem. Consequently, the results presented in Jl2| were obtained with this solver. 



4. Solving the Einstein/fluid system 

This section deals with some details of our numerical solution of the coupled 
Einstein/fluid equations, including the incorporation of source terms into our 
conservation laws, regularity and boundary conditions, and methods for calculating 
physical values for w in the ultrarelativistic regime. In addition, we describe the initial 
data and mesh structure we have used in our studies of critical phenomena in fluid 
collapse. Finally, we conclude the section with some remarks on how we have tested 
and validated our code. 



4-1. Time integration 



In section 2.4 the fluid equations of motion were written essentially in conservation 
form, except that a source term, i[>, had to be included. While this source term 
clearly breaks the strict conservation form of the equations, it can be self-consistently 
incorporated into our numerical scheme by using the method of lines to discretize 
space and time separately. Specifically, the discretized fluid equations become 
dq.i 1 



dt 



<r 2 X¥" 



i+l/2 



-1/2. 



(76) 
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where q 4 is the cellular average of q, F i±1 / 2 are the numerical fluxes defined by ([74]), 
and X = a/a, as previously. These equations can be integrated in time using standard 
techniques for ODEs. In particular, Shu and Osher |37| have investigated different 
ODE integration methods, and have found that the modified Euler method (or Huen's 
method) is the optimal second-order scheme consistent with the Courant condition 
required for a stabile evolution. We briefly digress to define this scheme for a general 
set of differential equations of the form 

(77) 

where L is a spatial differential operator. Let q™ be the discretized solution at time 
t = nAt, and L be the discretized differential operator. The modified Euler method 
is a predictor-corrector method, with predictor 

q*=q» + AiL(q"), (78) 

and corrector 

q ,l+1 = i(q" + q*) + ^AtL(q*). (79) 

Again, we note that At is subject to a Courant (CFL) condition, which can be deduced 
empirically or possibly from a linearized stability analysis. 

Particularly in comparison to the treatment of the fluid equations, numerical 
solution of the equations governing the geometric quantities a and a is straightforward. 
As discussed previously, the lapse, a, is fixed by the polar slicing condition (|l7]), while 
a can be found from either the Hamiltonian ( |l8| ) or momentum ( ft6| ) constraints. We 
have used discrete, second-order, versions of both equations for a, and have obtained 
satisfactory results in both cases (the polar slicing equation is likewise solved using 
a second-order scheme.) In general, however, (and particularly on vector machines) 
solution via the momentum constraint yields a far more efficient scheme, and we thus 
generally use the momentum equation to update a. 



Full details of our numerical scheme arc presented in Appendix B 



4-2. Regularity and boundary conditions 

In the polar-areal coordinate system, the lapse "collapses" exponentially near an 
apparent horizon, preventing the t = constant surfaces from intersecting the physical 
singularity which must develop interior to a black hole. As the slices "avoid" the 
singularity, elementary flatness holds at the origin for all times in the evolution, giving 

a(0,t) = l. (80) 

At each instant of time, the polar-slicing condition ([l7]) determines the lapse only 
up to an overall multiplicative constant, reflecting the reparameterization invariance, 
t — > t(t), of the polar slices. We chose to normalize the lapse function so that as r — > oo, 
coordinate time corresponds to proper time. On a finite computational domain, and 
provided no matter out-fluxes from the domain, this condition is approximated via 



aa 



1. (81) 



In spherical symmetry the fluid flows along radial lines, and given that there are 
no sources or sinks at the origin, we have that v(0, t) — 5(0, t) = 0. Thus 

n(0,t) = $(0,t) =r(0,t). (82) 
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Regularity at the origin further require r, IT and $ to have even expansions in r as 
r -> 0: 

T(r,t)=T (t)+r 2 r 2 (t) + O(r i ) (83) 

II(r, t) = n (t) + r 2 II 2 (£) + 0(r 4 ) = ro(t) + r 2 n 2 (<) + 0(r 4 ) (84) 

$(r, *) = $„(t) + rH 2 (t) + 0(r 4 ) = T (t) + r 2 $ 2 (<) + 0(r 4 ) (85) 

On our radial grid r^, i = 1, 2, • • • N, we use these expansions to compute grid-function 
values defined at r = r± = in terms of values defined at r — r 2 and r = r$. 
Specifically, once the values <I> 2 and $3 have been updated via the equations of motion, 
we compute $1 using a "quadratic fit" based on the expansion (85): 

*x = *^-*f a . (86) 

We then set ITi = $1. 

At the outer boundary we apply out-flow boundary conditions, which in our case 
are simply first-order extrapolations for II and 

$jv = $n-i njv=njv-i. (87) 

In addition, two ghost cells (r = r/v+ij r = r N+2) are added at the outer edge of the 
grid for ease in coding the cellular reconstruction algorithm [jl9| . These ghost cells are 
also updated with first-order extrapolation. 



4.3. Floor 

The fluid model is a continuum approximation, and, at least naively, the fluid equations 
become singular as p — > 0. In these evacuation regions, both the momentum and 
mass density are very small, and therefore the velocity — which loosely speaking is 
the quotient of the two — is prone to fractionally large numerical errors. These errors 
then often result in the computation of unphysical values for the fluid variables, such 
as supraluminal velocities, negative pressures or negative energies. (In addition, of 
course, our code must contend with the usual discretization and round-off errors 
common to any numerical solution of a set of PDEs.) At least from the point of view of 
Eulerian fluid dynamics, it seems fair to say that a completely satisfactory resolution 
of the evacuation problem does not exist. In the absence of a mathematically rigorous 
and physically acceptable procedure, we adopt the ad hoc approach of demanding that 
p > everywhere on the computational domain, i.e. we exclude the possibility that 
vacuum regions can form on the grid. In terms of our conservation variables q, this 
requirement becomes IT > and $ > 0. In a wide variety of situations, our numerical 
solutions of the fluid equations naturally satisfy these constraints. However, the critical 
solutions for "stiff" equations of state (r > 1.9) develop extremely relativistic velocities 
(W > 10 6 ) in regions where p is small [12], and we are unable to solve the fluid PDEs 
in these cases without imposing floor (or minimum) values on q. Specifically, at each 
step in the integration we require 

n > 6, <f>>6, (88) 

where the floor 8 is chosen to be several orders of magnitude smaller than the density 
associated with what we feel are the physically relevant features of the solution — a 
typical value is 8 = 10~ 10 . The floor is often applied in regions where IT and $ differ 
greatly in magnitude, and discretization errors can easily lead to the calculation of a 
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negative value for either function. For example, the floor may be applied to the "in- 
going" function in a region where the fluid is overwhelmingly "out-going." In these 
cases, the effect of the floor is dynamically unimportant. However, the floor may be 
invoked in other cases, where its effect on the dynamics is less certain. 

Given the ad hoc nature of this regularization procedure, the crucial question is 
whether the floor affects the computed solutions in a substantial way. We investigated 
this question by comparing critical solutions for r = 2 (the most extreme case) which 
were calculated with two distinct floor values: 5 = 1CP 8 , and our usual 5 = 10~ . 
The two solutions appeared identical, and identical mass-scaling exponents | ft2| were 
calculated. However, we note that the use of a floor makes estimates of the maximum 
Lorentz factor attained in the critical solutions unreliable because the largest velocities 
occur in regions where the floor is enforced. 



4-4 ■ Calculating the velocity 

The simple expression ( p6| ) for v in terms of q, when used naively with finite precision 
arithmetic, can result in the computation of unphysical, supraluminal velocities. For 
example, when searching for critical solutions we routinely calculate fluid flows with 
W > 10 3 . Thus, when calculating v from the quotient (^6|), small numerical errors 
can easily conspire to give |u| > 1, rather than the correct |u| > 0.999999. On the 
other hand, the combination 

X = W 2 v (89) 

is insensitive to small numerical errors, and provides a better avenue for calculating v 
from the conservation variables. From the definition ([l4|) of S we have 

The velocity can then be calculated from \ using 
1 



v= — [^l + Ax 2 -l) . (91) 
To the limit of machine precision, v is then in the physical range —1 < v < 1. When 



X *C 1, we calculate v from a Taylor expansion of (91), although ( |26| ) could also be 



used. We also use x when calculating w from q for the ideal-gas EOS 



4.5. Grid 

The black-hole-threshold critical solutions — which are our primary focus — are 
generically self-similar, and as such, require essentially unbounded dynamical range 
for accurate simulation. Thus some sort of adaptivity in the construction of the 
computational domain is crucial, and, indeed, the earliest studies of critical collapse jlj 
used Berger-Oliger adaptive mesh refinement |3£| to great advantage. However, in 
contrast to the early work, we know (at least schematically) the character of the critical 
solutions we seek, and thus we can, and have, used this information to construct a 
simple, yet effective, adaptive grid method. (Our approach is similar in spirit to that 
adopted by Garfinkle |3£j in his study of scalar field collapse.) Specifically, at any time 
during the integration our spatial grid has three distinct domains: the two regions near 
r = and r = r max have uniform grid spacings (but the spacing near r = is typically 
much smaller than that near the outer edge of the computational domain) , and the 
intermediate region has grid points distributed uniformly in log(r) (see figure |3|). As 
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Figure 3. Illustration of the re-meshing algorithm used in investigations of 
critical collapse. The grid spacing Ar is shown as a function of r on a log-log plot. 
The solid line represents the initial grid, the dotted line shows the grid spacing 
after the first addition of points near the origin, and the dashed line shows the grid 
spacing after the second regridding. Note that the grid spacings near the origin, 
and near the outer edge of the computational domain are uniform (horizontal 
lines). At each regridding cycle, the grid spacing near the origin is halved, and 
the new points are matched smoothly onto the previous grid. A critical evolution 
may involve more than 20 regriddings, although only a small number of points 
(50-150) may be added at a time. 

a near-critical solution propagates to smaller spatial scales, additional grid points are 
added in order to maintain some given number of grid points between r = and some 
identifiable feature of the critical solution. For example, we typically require that at 
least 300 or so grid points lie between the origin and the maximum of the profile of 
the metric function a. 

The primary advantage of this gridding scheme is that it is simple to implement, 
and yet allows us to resolve detail over many length scales: the ratio of the grid spacing 
at the outer edge to the spacing at the origin is typically 10 10 -10 13 at the end of an 
evolution. The primary disadvantage of this scheme is that it is specialized for critical 
collapse, and cannot be used for more general physical problems. 

4-6. Initial data for critical solutions 

We expect that the critical solutions in fluid collapse will be universal, in the sense that 
any family of initial data which generates families that "interpolate" between complete 
dispersal and black hole formation, should exhibit the same solution at the black hole 
threshold. We have thus focused attention on a specific form of initial data, which 
generates initially imploding (or imploding/exploding) shells of fluid. Specifically, the 
energy density in the shells has a Gaussian profile, 

r = r cxp [-(r - r ) 2 /A 2 ] + K, (92) 
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where the constant K — typically of magnitude 10 _6 r o — represents a constant 
"background". It should be note that this background is used only in setting the 
initial data, and is not held fixed during the evolution — in particular K is not a floor 



as discussed in section 4.3. The shells are either time-symmetric, or have an initial 
inward velocity which is proportional to r. Critical solutions were found by fixing r a 
and A, and then tuning the pulse amplitude t . 

4.7. Tests 

When developing a code such as the one described here, there are a number of tests 
which should be passed in order to provide confidence that the algorithm is producing 
reliable results. Perhaps most fundamental of these is the convergence test, which 
generally demonstrates that the numerical method is consistent and has been correctly 
implemented, but which also provides an intrinsic method for estimating the level of 
error in a given numerical solution. For our high-resolution shock-capturing scheme, 
a general rule-of-thumb is that the convergence should be (apparently) second order 
where the flow is smooth, and first order at discontinuities, where the effects of the 
slope limiter become important. In addition, we can also expect first order convergence 
near extrema of q, since at these points, the slope, s, changes sign, and the minmod 
limiter gives a piece-wise constant reconstruction for q. A convergence test where 
these effects are apparent is shown in figure |J 

After the numerical algorithm has been correctly implemented, one often 
compares results from the code to known closed-form solutions. In the early stages of 
code development, we tested the shock-capturing algorithm in this fashion by solving 
initial data for a shock tube, and comparing the results with the known solution of 
the Riemann problem. While the shock-tube provides a good test of the fluid solver, 
the test is done in Minkowski space with slab symmetry, and can probe neither the 
implementation of the geometric factors in the fluid equations, nor the discretized 
Einstein equations. A few general relativistic fluid systems can be solved exactly, and 
have traditionally been used to test new codes. These include static, spherical stars 
(Tolman-Oppenheimer-Volkoff), spherical dust collapse (Oppenheimer-Synder) and 
"cosmological" tests with a Robertson- Walker metric. In our companion paper |Q, 
we advocate the use of perfect fluid critical solutions as an additional test problem; 
one which involves both dynamic gravitational fields and highly relativistic fluid flows. 
Thus we consider the ultimate test of our full GR/fluid code to be the dynamical 
calculation of self-similar perfect fluid critical solutions, which can then be compared 
to solutions computed directly (but also numerically!) from a self-similar ansatz p3|. 

Appendix A. Characteristic structure 

In this appendix we calculate the Jacobian matrix A for the relativistic fluid equations, 
and then compute the associated eigenvalues and right eigenvectors. The flat-space 
components of A are 

1 BP 
All = -(l + 2v-v 2 )+(l-v 2 ) m (A.l) 

1 BP 
-Ai2= -2 {v + 1)2 + {1 ~ v2) d$ (A - 2) 
1 BP 

A2i^^(v~lf + (v 2 -l)— (A.3) 




Figure 4. Illustration of some of the convergence properties of the solution 
algorithm discussed in the text. Here we evolved a time-symmetric shell of fluid 
(r = 1.3) using uniform grids with three different resolutions: Ar = h, 2h and 4/i. 
Convergence is investigated by comparing the solutions obtained using the three 
distinct discretization scales. In frame (c), the solid line is [ji^ — t^) and the 
dotted line is 4 (t^ — T2h), where the subscript on r indicates the grid spacing for 
a particular solution. When the convergence is second order, the two lines should 
(roughly) coincide, while when the convergence is first order, the amplitude of 
the dotted line should be twice that of the solid line. As expected, we see that 
the convergence is not second order at the shock. (Of course the whole notion 
of convergence at a discontinuity fails, as the notion of Richardson expansion 
requires smooth functions.) However, we also can see that the convergence is only 
first order at the extrema of q — at these points, the slope changes sign, and the 
minmod limiter produces a first-order reconstruction. Frame (d) shows a more 
detailed view of a portion of the data displayed in (c). For context, we also show 
r in frame (a) and v in frame (b). 
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A 22 =U-l + 2v + v 2 ) + (v 2 - l)° (A.4) 



i(~l + 2. + . 2 ) + (^~l)§|, 

and the partial derivatives of P are easily found from (|3^). The eigenvalues A± of A 
are the two roots of the quadratic equation 

A 2 - (-An +-4 2 2)A + det^ = 0, (A.5) 

and the right eigenvectors are 

1 



V 



(A.6) 



where 



Y ± . ^Ai. (A.7) 

-^12 

If the eigenvalues are numerically degenerate owing to the limitations of finite precision 
arithmetic, we set A± = 0. When r = 2, the eigenvalues and eigenvectors become 
simply 

A±=±l, r+=(j), r_ = (M. (r = 2) (A.8) 



Appendix B. Implementation Details 



The origin in spherical symmetry requires additional care because powers of 1/r 
appear in the flux and source terms. One particular difficulty results from the 
partial cancellation of the source term, 2aP/(ar), with the pressure term in the 
flux. Numerically this cancellation is not exact, and this non-cancellation can induce 
large errors near the origin. We therefore modify the difference equations in order to 
eliminate the offending term. We first decompose the numerical flux into two parts 
f«andf( 2 ): 

P 

- P 



f(i) = 



|(n-$)(i + u) 
|(n- 



f(2) = 



so that f = f' 1 ) + f( 2 \ We then rewrite the conservation equations 
new fluxes as 



(B.l) 
with these 



where the new source term S is 

e 



s = 



-e 



The numerical flux function F is similarly decomposed: F = F^ 1 ) + F^ 2 ', with 
17(1) _ 1 



P (2) = 
i+1/2 



1 



fW(qf +1/2 ) + f (1) (q[ + i/ 2 ) - E l^l A ^r, 
f (2) (^ + i/ 2 ) + f (2) (qm/ 2 )]- 



(B.2) 
(B.3) 

(B.4) 
(B.5) 



The finite-differencing of the flux terms is adapted so that the derivatives have 
the correct leading order behavior near the origin. From the regularity conditions 
discussed in section 4.2 we have 

„3 



lim r 2 Af (1) 



lim Xf (2) cx constant, 



(B.6) 
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and we thus write the discretized equations of motion as 



dcu 
dt 



i+l/2 



f r 2 XF (l) 



i-1/2 



i+l/2 'i-1/2 
(^ F(2) ), + 1 /2-(^ F(2) 



i-1/2 



1*. (B.7) 

r i+l/2 — r i-l/2 

The geometric equations arc differenced using standard second-order finite- 



difference techniques. The momentum constraint is 
da t 2 

= /TTT-rv-r 

dt 



27rr i a i a? (11* - 



and is integrated using the modified Euler method described in section 4.1 
discretized polar slicing condition (|l7|) in discrete form is 



(lna)« +1 = (lna)? 



27rr((n- + 



2r 



1 



(B.8) 
The 



(B.9) 



i+l/2 



where all of the basic variables — a, II, &,v and P — in the {} braces are evaluated at 
r i+i/2 using an arithmetic average. 

Finally, the overall flow of an integration step is as follows: 

(i) Begin with the data for time t = t n : {IP, $™, P™, v n , a n , a 71 }. 

(ii) Reconstruct cells for {q £ ,q r }, and calculate the numerical fluxes F(q £ ,q r ). 

(iii) Perform the predictor step of the modified Euler method, obtaining {II*, $*, a*}, 
then calculate {P*,v*}, and integrate the slicing condition to determine a*. 

(iv) Reconstruct cells for {q^*,q r *}, and calculate the numerical fluxes F(q f *,q r *). 

(v) Perform the corrector step of the modified Euler method, obtaining 

jra+l (j>«+l a n + l 

condition to determine a n+1 

(vi) Check the regridding criteria, and regrid if necessary. 



{n"+ 1 ,$ ,l+1 ,a™+ 1 }, then calculate {P n+1 ,v n+1 }, and integrate the slicing 
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